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Abstract. Single molecule scanning tunneling spectroscopy (STS), with dephasing 
due to elastic and inelastic scattering, is of some current interest. Motivated by this, 
we report an extended Huckel theory (EHT) based mean-field Non-equilibrium Green's 
function (NEGF) transport model with electron-phonon scattering treated within the 
self-consistent Born approximation (SCB A) . Furthermore, a procedure based on EHT 
basis set modification is described. We use this model to study the effect of the 
temperature dependent dephasing, due to low lying modes in far-infrared range for 
which TujJ <C ksT, on the resonant conduction through highest occupied molecular 
orbital (HOMO) level of a phenyl dithiol molecule sandwiched between two fec- 
Au(lll) contacts. Furthermore, we propose to include dephasing in room temperature 
molecular resonant conduction calculations. 
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Figure 1. (color online) Schematic diagram. Scanning Tunneling Spectroscopy (STS) 
setup for probing conduction through phenyl/benzene dithiol (PDT/BDT) molecule 
bonded to Au(lll) substrate. Voltage is applied at the tungsten W(110) tip. (a) 
A schematic showing the self-consistent procedure of solving EHT-NEGF equations 
for coherent transport [22] • (b) A schematic detailing the self-consistent procedure of 
solving EHT-NEGF-SCBA equations for quantum transport with elastic and inelastic 
dcphasing. 



1. Introduction 

We present an extended Hiickel theory (EHT) based self-consistent mean-field transport 
model for incoherent single molecule spectroscopy. Although simple, this approach 
is worth-pursuing because EHT is numerically inexpensive and overcomes many of 
the short-comings of the more sophisticated theories [1] and has been used for 
nanoelectronics [2J. Furthermore, the simplicity and the exponential dependence of the 
Slater type orbitals basis set allows an intuitive solution to tip modeling as described in 
this paper. The proper tip modeling is important because if condensed matter like short 
ranged orbitals are used, one has to bring the tip very close to the molecule [about 4-5A] 
to obtain the experimentally observed current levels. This artificially results in higher 
than the actual voltage drop across the molecule resulting in not only wrong Laplace's 
potential but Hartree potential is also overestimated. 

Furthermore, the electron-phonon scattering in molecular structures on different 
metallic substrates has been studied quite extensively at low temperature experimentally 
US m El E] and theoretically [3 El U UHl Ell CL1 13] in the form of off-resonant inelastic 
tunneling spectroscopy (IETS). In IETS, the molecule acts as the scattering region and 
the distinct peaks in the second derivative of I-V characteristics identify the molecular 
vibrations. Moreover, in electron-phonon scattering, not only electron exchanges energy 
with the phonon bath but also its phase is irreversibly lost due to phonon degree of 
freedom. If the phonon energy (hoS) is less than a characteristic energy (e.g. ksT), 
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this inelastic dephasing process may be approximated by an elastic dephasing event for 
numerical convenience. Since all the energy channels are coupled in inelastic dephasing, 
solving for it is computationally challenging. 

Apart from this, the phonon modes of a molecule bonded to a substrate are 
different from gas phase. Organic molecules in gas phase have vibrational energies in the 
infrared range, i.e. 100-350meV. However, experiments using IETS [5j, high resolution 
electron energy loss spectroscopy (HREELS) [13] and low temperature transport [H] 
suggest that bonding to a substrate leads to low lying modes in the far-infrared range. 
Dominant modes having vibrational energies less than about lOmeV have been reported 
theoretically [HIE]. These far-infrared modes are the result of the combined motion of the 
molecule with respect to the contacts. Most of these infra-red modes are longitudinal 
modes and couple well with carriers 0, |9j. Hence, the low lying modes, for which 
hu <C ksT at room temperature, are the most important phonon modes for molecules 
bonded to substrates. In a previous study [15J, we report transport calculations at room 
temperature for a styrene chain bonded to n ++ -H:Si(001) substrate. In Ref. [15], we 
propose to include dephasing due to low lying modes in order to theoretically reproduce 
experimental results. 

Furthermore, in context of molecular spectroscopy, there are four sources of 
broadening: (1) contact broadening which is about 0.1-0.3eV for a good contact (2) 
the Fermi function broadening of the two contacts (3) broadening due to Hartree 
self-consistency which is usually small in STS setup due to tip being about lnm 
away from substrate and (4) broadening due to electron-phonon scattering which is 
temperature dependent due to phonon population given by Bose-Einstein factor as 
l/[exp(huj/kBT) — 1]. The main focus of this work is the broadening due to this fourth 
source in transport through a phenyl/benzene dithiol (PDT/BDT) molecule bonded to 
Au(lll) substrate as shown in Fig. 1. The tip configuration used is tungsten W(110) 
having a workfunction of about 5.2eV. 

The paper is divided into five sections. In Sec. II, we describe the mean-field 
self-consistent EHT-NEGF-SCBA model for studying transport with dephasing due to 
electron-phonon scattering. Furthermore, details are provided about the electrostatic 
calculation. Tip modeling is an important aspect of the calculations and is discussed 
in Sec. III. In Sec IV, we discuss the results and report the temperature dependence 
of the transport quantities on the electron-phonon scattering due to low lying modes. 
Finally, we provide the conclusions in Sec. V. The molecular visualization is done using 
Gauss View [16J. 

2. Formalism 

We use the single particle non-equilibrium Green's function formalism (NEGF) in the 
mean field approximation using non-orthogonal basis to model the quantum transport. 
We define the time-retarted non-equilibrium Green's function in non-orthogonal basis 
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Figure 2. (color online) Block diagram for transport with elastic dephasing. 
For elastic dephasing, all the energy channels are independent. Therefore, energy 
dependent quantities can be calculated independently for each energy point with the 
convergence criteria shown. Equilibrium and non-equilibrium quantities for each bias 
point are calculated as shown. The green and blue blocks are color coded to match with 
the blocks in Fig. l(a,b). For inelastic dephasing, since energy channels are coupled, 
one has to use convergence criteria on the calculated energy dependent quantities over 
the whole energy range. 



as PUCE]: 

G(E, V tip ) = [(E + iO + )S -H- E fc S - S lj2 - E^ 1 (1) 

where S is the overlap matrix and H = Hd + Ud- Hd is the device Hamiltonian for an 
isolated molecule expressed in EHT [19]. S 12 are the self-energies of W(110) tip and 
Au(lll) contact respectively - calculated recursively using EHT as well. Since, there 
is no vacuum reference in EHT parametrization, we use input from experiments or 
other sophisticated theoretical methods to set the reference. Ef c is used as a parameter 
to achieve this objective. It is taken 2eV in this paper to match the experimentally 
observed peak in conductance spectrum. Ud incorporates both Laplace's potential due 
to the applied tip voltage (Vnp) and Hartree potential due to the non-equilibrium density 
of electrons in the molecular channel. The contact self-energies S 12 are related to the 
broadening functions as Ti j2 = 2(Si j2 — The tip self-energy (Ei) is further discussed 
in Sec. III. The Au(lll) self-energy (S 2 ) has been discussed in Ref. [22]. The spectral 
function is defined as A = i(G — G^), which is related to the density of states of the 
molecular region as D(E) = 2(for spin) x tr(AS) /2tt . The electron correlation function 
G n (= — zG < ) is given as: 

where £q n 2 (= — ^^2) are the contact inflow functions defined as r 12 / li2 and S* n (= 
—iLf) is the scattering inflow function, f 12(E) are the contact Fermi functions given 
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as 1/[1 — exp((E — yU 12 )//c B T)]. // li2 are the chemical potentials of the two contacts given 
as /ii = fi a — eVup and fi 2 = /^o respectively - /i D is the equilibrium chemical potential. In 
this model, contacts are assumed to remain in equilibrium. This is a valid assumption 
for metallic contacts and degenerately doped semi-conducting contacts. However for 
the moderately doped or lightly doped semi-conducting contacts, band bending due to 
V t i P makes these calculations computationally expensive. The hole correlation function 
is defined as G p = A — G n . The electron density matrix is then given as: 

p = _L r dE G n (E) (3) 

Z7T J-oo 

The total number of electrons is computed as N = 2(/or spin) x tr(pS). Finally, the 
I-V characteristics are computed as follows [171 120] : 

Ii(V Up ) = 2(/or spin) x - / dE tr(Y!?A - T t G n ) (4) 

tl J — oo 

Since scattering self-energy (S s ), G, G n and S* n depend on each other, we solve these 
four quantities self-consistently along with the Hartree self-consistent loop as shown in 
Fig. 1(b). The flow diagram for the self-consistent procedure at each bias is shown in 
Fig. 2. For SCBA, it can be shown that the current through the scattering contact is 
always zero. 

The scattering inflow function outflow function and broadening 

function (T s ) are defined as [21]: 

Y™(E) = I" ^^D em (hu)SG n (E + hu)S + D ab (huj)SG n \E - huj)$h) 
Jo 2n 

T^\E) = r ^lD em {huj)SG p {E - hu)S + D a \huo)SG p {E + huo)$>) 
Jo 2ir 

F S (E) = Ef + (7) 

= r^^D em (huj)S[G n (E + huj + G p (E-huj)}S (8) 
Jo 2i 

+D ab (hu)S[G n {E - huj) + G P {E + hu)]S (9) 

where D em (huj) = (N + l)D (fiu) and D ab (hu) = N D (hu>) are the emission and 
absorption dephasing functions respectively. N is the equilibrium number of phonons 
given by Bose-Einstein statistics as l/[exp(hu /ksT) — 1]. Eq. 5 implies that electron 
inflow at a particular energy E is dependent on the electron correlation function G n at 
energy E + hu for emission and E ~hu for absorption. This relationship is reversed for 
the electron outflow as in Eq. 6. The overall broadening due to this scattering process 
is given by T S (E). Furthermore, D (hu) is a fourth ranked tensor and is related to the 
electron phonon interaction potential (U) as < i,j\WU\k,l >. Thus, D (hu) can be 
different for different molecular orbitals and is explicitly written as D Q (i, j; k, I; huj). Eqs. 
5-7 are modified from those reported in Ref. [21] for non-orthogonal basis as discussed 
in Ref. HE]. 
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High energy phonon limit: For hu 3> kgT, N pa and N + 1 ~ 1, hence D em pa D 
and D ab pa 0. Thus, Eq. 7 becomes: 

TJE) = f°° ^lD (hu)S[G n (E + hu) + G P (E - hu)]S (10) 

JO 27T 

Thus, contribution due to absorption is negligible in this case because not enough 
phonons are available to be absorbed. Moreover, since T S (E) is temperature 
independent, the broadening due to electron-phonon scattering for hu ^> ksT is 
approximately temperature independent. 

Low energy phonon limit: For hu <C ksT, by keeping the first term in Taylor series 
expansion of exp(hu ksT) pa 1 + huf/ksT, one obtains iV + 1 pa N pa ksT/hu. This 
leads to, 

d(hu) 



and 



FJE) = I" ^l[D em (hu) + D ab (hu)]SA(E)S (11) 



1 rp 

D em (hu) pa D a \hu) = N D (hu) Pa J?—D (hu) (12) 

hu 



We finally obtain T S (E) as, 

r«2yyM (13) 

Jo 2ti hu 

D 

where D is referred to as dephasing strength. As shown above, the broadening due to 
elastic dephasing is directly proportional to temperature. 

Furthermore, the real part of T, S (E) is computed by taking a Hilbert transform of 
the imaginary part as below, 

E S (E) = - r dy ~ Ts{y)/2 - iT s {E)/2 (14) 
7T J-oo E — y 

Calculating the real part of Yi s (E) is computationally expensive and leads to convergence 
issues. Previously, it has been reported [10] that this real part contributes on the order 
of meV. We therefore ignore it in this paper. We phenomenologically approximate the 
dephasing strength (D) by a constant. This approximation treats T\ and T 2 time in 
an average manner and seems to reproduce experimental results [15]. Treating D as a 
constant is similar to the dephasing strength used in a Biittiker probe. However, the 
mathematical framework of SCBA is of course different from that of Biittiker probe 
method. Furthermore, a rigorous calculation of dephasing functions and dephasing 
strength can be readily incorporated in our model. Apart from this, SCBA has some 
limitations, e.g. it does not take care of multiple phonon processes implicitly for a single 
scattering event and it does not include crossed diagrams which correspond to successive 
emission and re-absorption of phonons by the same electron. 

We use a finite-element solver [23] for solving Laplace's equation in three dimension 
to obtain potential drop across molecule due to Vu p . This method also provides the zero 
bias band line up potential due to the Fermi level mismatch of O.leV between Au(lll) 
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Figure 3. (color online) Tip modeling. The calculated current (I) for Vu p = XV 
as a function of tip-molecule distance (z) showing exponential dependence of current 
on z. The extracted barrier height from this I-z plot is 4.3eV, whereas with original 
parameters it is approximately lOeV. The calculated apparent barrier height (<f> ap ) is 
close to the work functions of the contacts being used. Furthermore for calculating 
Laplace and image potential due to Vu p , a 20nm diameter tip is used as shown in the 
inset. The tip is 7A away from the molecule in this visualization. 



and W(110) tip, whose work functions are 5.1eV and 5.2eV respectively. We take proper 
tip shape into account. A 20nm tip diameter is assumed as show in the inset of Fig. 
3. The Hartree potential for the molecule is approximated via the complete neglect 
of differential overlap method [221 1211 • Image effects are incorporated to ensure that 
self-consistency is not over-estimated. 

3. Tip Modeling 

For a good tip which can give atomic resolution, the last atom at tip apex dominates 
STS. Therefore, we calculate the overlap between molecule and tip for this single 
tungsten atom. This assumption is for a "good" tip without any adsorbate or multiple 
tip structures which may give ghost images in scanning tunneling microscopy (STM). 
We use W(110) configuration since it is the lowest energy state. Furthermore for working 
tip heights, electronic effects due to the wave-function overlap between the tip atoms 
and the molecule become important. This overlap gives exponential dependence of the 
decay of tunnel current (I) with tip height (z). Within WKB approximation, apparent 
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Figure 4. (color online) Incoherent resonant conduction through the highest occupied 
molecular orbital (HOMO) level with elastic dephasing due to low energy electron- 
phonon scattering. The dephasing strength (D) is directly proportional to the 
temperature (T). Thus, the additional broadening due to elastic dephasing decreases 
with decreasing temperature in a linear fashion, (a) I-V characteristics showing 
broadening of spectroscopic features due to dephasing. (b) dl/dV-V characteristics 
providing an alternate and more detailed view of broadening of spectroscopic features, 
(c) Full width half maximum (FWHM) of the first conductance peak corresponding to 
the HOMO level showing linear dependence on temperature. The dephasing strength 
(D) is 0.15eV 2 at 300K to give experimentally observed broadening. 



barrier height (<f> ap ) can be calculated from the slope of the log(T)-z plot. This apparent 
barrier height is usually found to be close to the workfunctions of the materials being 
used - Au(lll) in our case. We find that the apparent barrier height is about lOeV with 
unmodified EHT parameters, which is unphysical. We report a scheme for modifying 
basis set of tip atom to get correct apparent barrier height. One does not need to modify 
the basis set of atoms in the molecule because they are already optimized for gas phase. 
Similar approach has been used before in the context of scanning tunneling microscopy 
[25J. With the modified EHT parameters of the s-orbital basis set [26] for tungsten 
atom, we obtain <ft ap = 4.3eV from the calculated I-z plot with Vu p = IV as shown in 
Fig. 3. 

The above parameter modification affects Si and Hi, provided the tip self-energy 
Si is defined as [(E + iO + )Si — Hi]g s i((E + iO + )S{ — H\), where g s i is the surface Green's 
function for the tip. We assume a constant g sl , calculated from the density of states 
at equilibrium chemical potential of bulk tungsten Dw(fJ-o), z.e. g s i = —inDw(tio) - a 
commonly used approximation [see e.g. Ref. [15] and references there-in]. However, 
Si is still energy dependent and this energy dependence should in fact capture all the 
barrier tunneling effects. However, it would not be able to capture the features due to 
peculiar electronic structure of tip as in Ref. [2~T] . 
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4. Results 

In this section, we apply the proposed EHT-NEGF-SCBA model to incoherent transport 
with elastic dephasing through the highest occupied molecular orbital (HOMO) level of 
a PDT/BDT molecule bonded to a Au(lll) substrate through a single Au-S bond and 
probed using a tungsten W(110) tip as shown in Fig. 1. Au-S bond length is the 
standard 2.53A and tip molecule distance is KM. The dephasing strength (D) at 300K 
is taken as 0.15eV 2 to obtain experimentally observed broadening. Dephasing strength 
at any other temperature T is then calculated following Eq. 11. 

The calculated I-V characteristics are shown in Fig. 4(a). For coherent conduction 
through PDT/BDT molecule at 300K (D=0), the I-V curve has a sharp rise. Fig. 4(b) 
provides an alternate view showing the dl/dV-V characteristics. The full width half 
maximum (FWHM) of the conductance peak is about 0.2eV. This broadening is due 
to (1) the finite life-time of electrons in the molecule as a result of the coupling with 
the contacts given by the broadening functions Ti^, (2) contact Fermi functions giving 
about 5.4/c#T broadening, and (3) Hartree potential, which is small because tip is 10 A 
from molecule. 

With D = 0.15eV 2 at 300K, the I-V curve gets smeared out due to spectral 
broadening induced by elastic dephasing and overall current level decreases as shown in 
Fig. 4(a). Furthermore, current starts increasing earlier and hence is slightly higher than 
D=0 case for V tiv < 1.45V. Correspondingly, the dl/dV plot [see Fig. 4(b)] has a higher 
conductance value for 1.45V > Vu p > 1.8V, (i.e. before and after the conductance 
peak). The FWHM of this conductance peak is about 0.8eV as shown in Fig. 4(c). The 
additional 0.6eV broadening is due to elastic dephasing. Such a high broadening has 
been observed experimentally in Ref . [28J . Apart from this, the comparison of I-V and 
dl/dV-V plots for D=0 and D = 0.15eV 2 at 300K raises two important questions (1) 
why the overall current value decreases and (2) how does the temperature dependence 
of dephasing strength affect the FWHM of conductance peaks. 

The answer to the first question depends on the energy dependence of the contact 
broadening functions r 12 . After including dephasing, the spectral weight increases in 
the energy range where it was smaller before and vice versa. If the energy region, where 
the spectral weight is transferred, has larger broadening functions r^, it results in a 
higher current. We analyze the energy dependence of r 12 and find that r 12 decrease as 
a function of decreasing energy, which corresponds to positive Vu p . Thus, the broadened 
molecular spectral function gives rise to higher current values for lower Vu p , where the 
broadening functions are larger. Similarly, current values are smaller for larger 
Vup, where r^a are smaller. Hence, the overall current level decreases with dephasing. 
However, in a previous study [15] with n ++ — H : 5i(001) — (2x1) contact, we conclude 
that the overall current level increases with dephasing. In Ref. [T5j, the HOMO level 
is near valence band-edge, where the contact density of states (DOS) is small resulting 
in a small current. After including dephasing, the HOMO level gets broadened and a 
portion of it is in energy range where the DOS is large. Thus, dephasing in this case 
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Figure 5. (color online) Effect of inelastic scattering on the off-resonant conduction 
at 4K assuming one phonon mode with D em = O.leV 2 and hui = IhmeV . (a) IETS 
signature is not noticeable in the I-V characteristics, (b) Inelastic scattering results 
in a step in dl/dV-V characteristics, (c) IETS spectra (d 2 I /dV 2 -V characteristics) 
showing a peak at 75meV corresponding to the phonon mode. Since Tiw ksT, 
D ab = D em e -hu;/kBT _ q ThuSj the IET s peak is due to em i ss i on f phonons and 

processes due to absorption of phonons are weak. 



results in an order of magnitude higher current. 

With decreasing temperature and dephasing strength, e.g. at 200K and 100K with 
D = O.leV 2 and D = O.OSeV^ 2 respectively, the I-V plots follow same trend as shown 
in Fig. 4(a,b). The calculated FWHM values of these conductance peaks are shown in 
Fig. 4(c) to quantify the effect of dephasing. FWHM is directly proportional to the 
temperature and dephasing strength. The dependence of FWHM of conductance peak 
on temperature and dephasing strength could become complex if temperature dependent 
structural changes are present. 

In Fig. 5, we report the inelastic tunneling spectroscopy (IETS) of PDT/BDT 
molecule at 4K. The purpose of these standard and well-established results is to convey 
the message that the model is able to handle calculations for inelastic dephasing as 
well. We assume a single phonon mode with energy hu = 75meV and emission 
dephasing function D em = O.leV 2 . The absorption dephasing function is then given 
as D ab = D em e ~ hu} / k BT _ pjg_ e^ a ~) gho^yg the j_y characteristics with no feature of IETS 
signal. Fig. 5(b) shows the dl/dV-V characteristics with a step in the conductance 
curve due to emission of phonons at eVu p = hu = 75meV. The same signal appears as 
a peak in the d 2 I/dV 2 spectrum as shown in Fig. 5(c). Furthermore, we include the 
Hartree self-consistency for elastic dephasing only. Since, we are studying the inelastic 
dephasing in the off-resonance regime, self-consistency is not important. 
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5. Conclusions 

We have reported an EHT-NEGF-SCBA model and have used it to study elastic 
and inelastic dephasing. It has been previously reported that low lying modes are 
important in molecular conduction at low temperature [HI [HI El [9] . We suggest that they 
should be important at room temperatures as well and should be included in transport 
calculations. Furthermore for huo <C ksT, inelastic dephasing can be approximated by 
elastic dephasing. Within this approximation, it can be shown that the broadening 
is proportional to temperature. Moreover, for Au contacts, the overall current level 
tends to reduce due to this dephasing. This could help bridge the gap between 
theory and experiment for resonant conduction along with other proposals [30l [31], 
where theoretically calculated currents are higher as compared to the experimental 
observations [29]. Finally, we present results for IETS to show that within the model, 
inelastic dephasing can be handled as well. 

It is a pleasure to acknowledge useful discussions with S. Datta and E. C. Kan. 
We thank F. Zahid and T. Raza for Huckel-IV 3.0 [22] codes and later for reviewing 
the manuscript. Computational facilities were provided by the NSF Network for 
Computational Nanotechnology. 
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